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1.  Introduction 


Dissipative  Particle  Dynamics  (DPD)  is  an  attractive  method  for  simulating  soft  condensed 
matter,  including  polymers,  surfactants,  and  colloids  at  the  mesoscale.  DPD  is  a  particle-based 
mesoscopic  method  that  can  be  understood  as  a  coarse-graining  technique  that  correctly  predicts 
the  hydrodynamic  nature  of  a  fluid.  The  general  idea  behind  DPD  is  that  particles  interact  with 
each  other  through  a  set  of  stochastic  differential  equations  with  conservative,  dissipative,  and 
random  forces. 

While  the  original  DPD  formalism  is  isothermal  ( 1 ,  2),  it  has  been  extended  to  the  isoenergetic 
(DPD-E)  (3,  4)  case,  which  is  of  practical  interest  for  simulating  materials  at  nonisothermal 
conditions.  One  challenge  requiring  special  consideration  is  the  numerical  integration  of  the 
stochastic  equations-of-motion  (EOMs).  The  most  commonly  applied  numerical  integration 
scheme  for  DPD  applications  originates  from  the  velocity- Verlet  (VV)  integrator  used  for 
molecular  dynamics  applications  and  is  extended  to  DPD  by  including  the  dissipative  and 
random  forces  in  the  overall  force  expression  (5).  This  approach  is  currently  used  in  the  Large- 
scale  Atomic/Molecular  Massively  Parallel  Simulator  (LAMMPS)  software  package  ( 6)  but  is 
limited  only  to  the  constant-temperature  DPD  method.  An  alternative  integration  scheme  using 
the  Shardlow-splitting  algorithm  (7)  separates  the  EOMs  into  stochastic  and  deterministic 
integration  steps.  It  is  readily  extendable  to  all  DPD  variants  and  has  been  found  to  be  the  most 
effective  integration  scheme  available  to  date,  especially  when  considering  DPD  simulations 
under  isoenergetic  conditions  to  attain  the  length  and  time  scales  necessary  to  model  complex 
systems. 

In  this  work,  we  present  a  general  framework  for  implementing  the  DPD-E  method  into  the 
highly  scalable  LAMMPS  simulation  software  to  efficiently  model  systems  under  isoenergetic 
conditions.  We  extend  the  current  LAMMPS  VV  integration  scheme  for  isothermal  DPD 
simulations  to  the  isoenergetic  case.  In  addition,  we  describe  the  implementation  of  the 
Shardlow-splitting  algorithm  (SSA)  to  enable  longer  time  steps  with  comparable  accuracy. 
Finally,  a  description  of  example  benchmark  problems  is  provided,  along  with  a  discussion  about 
the  tradeoffs  between  the  DPD  version  of  the  VV  and  the  SSA  integration  schemes  in  terms  of 
performance  and  stability. 
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2.  Methods 


2.1  Constant  Energy  Dissipative  Particle  Dynamics  (DPD-E) 

In  the  DPD-E  method  (3,  4)  particles  are  defined  by  a  mass  mi ,  position  r, ,  momentum  p, ,  and 
particle  internal  energy  Uj.  The  variation  of  d ut  is  taken  as  the  sum  of  the  two  contributions  that 
correspond  to  the  mechanical  work  done  on  the  system  d u™ech  and  the  heat  conduction  between 
particles  d uc°nd  .  For  constant-energy  and  constant-volume  conditions,  the  evolution  of  DPD 
particles  in  time  t  is  governed  by  the  following  EOM: 


where  r(/.  =  r  -  r  is  the  separation  vector  between  particle  i  and  particle  j ,  and  r  =  |r. 1 . 

P  p  . 

v,  =  — - —  is  the  relative  velocity  between  the  particle  pair.  yu  and  cr  are  the  friction 

v  mi  m 

coefficient  and  noise  amplitude  between  particle  i  and  particle  j ,  respectively,  d If  =  d  IE  are 

the  increments  of  the  Wiener  processes  associated  with  momentum  variations.  The  weight 
functions  a>D(r)  and  coR{r )  vanish  for  r>rc,  where  r  is  the  cutoff  radius,  and  are  typically 

chosen  as 

n)D(r)  =  [( ufi(r)]2  =  j(^  rc)  '  r  <  rc  (2) 

(  0  ,  r  >rc 

Similarly,  the  weight  functions  ojDg(r)  and  coRq(r)  are  defined  in  an  equivalent  manner  as 
equation  2.  0i  is  the  particle  internal  temperature  and  is  related  to  ut  through  a  mesoparticle 
equation  of  state.  Ktj  and  atj  are  the  mesoscopic  thermal  conductivity  and  the  noise  amplitude 
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between  particle  i  and  particle  j ,  respectively,  and  d  Wq  =  -d  Wq  are  the  increments  of  the 
Wiener  processes  associated  with  thennal  conduction.  The  conservative  force  Fl  is  given  as  the 
negative  derivative  of  a  coarse-grain  potential,  ufj ,  and  is  expressed  in  LAMMPS  as 


dutf 

drij 


AijO)R 


(3) 


where  Ay  is  the  force  constant.  Bonet  Avalos  and  Mackie  (3)  demonstrated  that  thermodynamic 
consistency  requires  the  following  fluctuation-dissipation  relations  to  be  satisfied 


coD(r)=  [o/(r )]" 
a]  =  2  kBKj 

(0Dq{r)=[(0Rq{r)] 


where  the  relevant  temperature  is  ©3  =  ^ 


x0*+0jj 


.  In  the  following  outline,  the  EOMs  are 


formulated  for  the  VV  integration  scheme  under  isoenergetic  conditions. 
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2.  Force  Calculation:  {F;  ,  where  F(  =  ^Fl7  +  F^5  +  F* 


j*i 


3.  For  i  =  1,. ..,7V 
a.  p'<-p, 
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Here,  the  EOMs  are  formulated  for  the  VV  integration  scheme  using  the  Shardlow-splitting 
algorithm  (VV-SSA)  (5). 


1 .  Stochastic  Integration :  For  all  i  -  j  pairs  of  particles 
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2.  Deterministic  Integration  #1:  For  i  = 

A  t 

a  p  ^p,+  — F, 
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2.2  Parallel  Implementation  of  DPD-E  Into  the  LAMMPS  Framework 

A  detailed  description  of  the  theoretical  foundations  and  parallelization  of  the  DPD-E  method 
using  the  VV  and  VV-SSA  integration  schemes  can  be  found  elsewhere  ( 8 ,  9).  In  this  report,  the 
relevant  source  code  modifications  that  were  required  to  implement  the  VV  and  VV-SSA 
integration  schemes  for  DPD-E  are  presented.  A  summary  of  the  code  modifications  is  provided 
in  appendices  A-G.  The  source  code  is  current  with  the  27  January  2014  version  of  LAMMPS 
and  is  available  upon  written  request  to  the  authors. 

To  implement  DPD-E  into  LAMMPS,  a  new  atom  style  was  created  to  handle  the  DPD  particle 
attributes  required  for  DPD-E  simulations.  In  addition,  new  pair  styles  were  created  to  compute 
the  DPD  forces,  and  new  fixes  were  created  to  integrate  the  DPD-E  equations  of  motion  through 
a  VV  and  VV-SSA  integration  scheme.  Compute  functions  were  created  to  monitor  the  DPD 
particle  attributes  as  a  simulation  progresses.  Each  new  feature  to  LAMMPS  is  described  in 
detail  in  the  following  sections. 

2.2.1  Implementation  of  the  DPD  Atom  Style 

For  all  isoenergetic  DPD  calculations,  a  new  atom  style  dpd  is  required  to  compute  and 
communicate  between  processors  the  DPD  particle  attributes  and  per- atom  arrays.  DPD  particles 
are  specified  in  the  LAMMPS  input  file  via  the  atom  style  command: 

atom_style  dpd 


Selection  of  the  dpd  atom  style  requires  the  DPD  internal  temperature  to  be  specified  in  the 
corresponding  LAMMPS  data  file  according  to  the  following  fonnat: 

•  Column  1 :  particle  id 
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•  Column  2:  particle  type 

•  Column  3:  Internal  temperature  (00  of  particle  i 

•  Columns  4-6:  Cartesian  coordinates  (x,  y,  z)  of  particle  i 

The  implementation  of  the  dpd  atom  style  to  LAMMPS  requires  a  new  AtomVecDPD  class  to  be 
created  in  order  to  compute  and  communicate  the  DPD  particle  internal  temperature  (dpdTheta), 
equation  of  state  flag  (eos),  conductive  energy  (uCond),  mechanical  energy  (uMech),  as  well  as 
the  differences  in  the  DPD  particle’s  conductive  energy  (duCond)  and  mechanical  energy 
(duMech)  between  two  subsequent  integration  steps.  The  AtomVecDPD  atom  style  class  derives 
from  the  AtomVec  parent  class  and  is  similar  to  and  modeled  after  the  existing  AtomVecAtomic 
(atom  style  atomic)  class. 

Implementation  into  LAMMPS  requires  modification  of  the  existing  atom  class  (LAMMPS 
source  files:  atom.h  and  atom.cpp)  and  creation  of  the  AtomVecDPD  class  (LAMMPS  source 
files:  atom_yec_dpd.h  and  atom_vec_dpd.cpp).  The  atom  class  code  is  modified  to  define  the 
necessary  pointers  to  the  DPD  attributes.  The  differences  between  the  modified  and  the  native 
LAMMPS  codes  are  shown  in  appendix  A.  Next,  the  new  AtomVecDPD  child  class  was  created 
by  copying  the  existing  AtomVecAtomic  class  (LAMMPS  source  files:  atom_vec_atomicdi  and 
atom_vec_atomic.cpp),  then  adjusting  the  data  file  format,  the  data  structure  sizes,  and  the 
communication  buffers  to  handle  the  additional  per-particle  arrays.  A  summary  of  the  differences 
between  the  AtomVecDPD  and  AtomVecAtomic  classes  is  given  in  appendix  B. 

2.2.2  Implementation  of  the  DPD  and  DPD/Atom  Compute  Commands 

For  all  isoenergetic  DPD  calculations,  it  is  useful  to  monitor  all  internal  properties  of  the 
particles  on  both  a  system  and  per-particle  basis.  The  ComputeDPD  and  ComputeDPDatom 
classes  were  created  to  monitor  these  properties  over  the  course  of  the  DPD-E  simulations. 

The  ComputeDPD  class  (LAMMPS  source  files:  compute_dpd.h  and  compute _dpd.cpp ) 
computes  the  total  particle  internal  conductive  and  mechanical  energies  by  summing  the  per- 
particle  energies.  In  addition,  the  particle  internal  temperature  of  the  system  is  computed  through 
a  harmonic  average  of  the  per-particle  internal  temperatures,  defined  as: 

N 

.  1  ST1  1 

dav9  =  M^J  (5) 

i= i 

where  N  is  the  number  of  particles  in  the  system.  The  ComputeDPD  class  is  accessed  through  the 
LAMMPS  input  file  via  the  compute  command: 

compute  1  all  dpd 

and  returns  a  vector  of  size  5  that  contains  the  following  particle  internal  properties: 
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•  The  total  conductive  energy  of  the  system 

•  The  total  mechanical  energy  of  the  system 

•  The  sum  of  the  conductive  and  mechanical  energy  of  the  system 

•  The  harmonic  averaged  internal  temperature  of  the  system 

•  The  number  of  particles  in  the  group 

The  ComputeDpdAtom  class  (LAMMPS  source  files:  compute _dpd_atom.h  and 
compute_dpd_atom.cpp)  accesses  the  per-particle  internal  energies  and  internal  temperature.  The 
compute  is  specified  through  the  LAMMPS  input  file  via  the  compute  command 

compute  1  all  dpd/atom 

and  enables  access  to  the  following  particle  properties: 

•  The  per-particle  conductive  energy 

•  The  per-particle  mechanical  energy 

•  The  per-particle  internal  temperature 

2.2.3  Implementation  of  DPD-E  With  the  VV  Integration  Scheme 

Implementation  of  the  isoenergetic  DPD  using  the  VV  numerical  integration  scheme  is  similar  to 
the  existing  isothennal  DPD  LAMMPS  implementation.  The  major  difference  is  the  calculation 
of  the  particle  internal  energies  within  the  dpde  pair  style  compute  command  and  the  integration 
of  the  internal  energies  within  the  fix  dpde  command.  The  new  pair  style  and  fix  commands  are 
presented  together  since  the  classes  interact  with  one  another  to  perfonn  all  stages  of  the  VV 
integration. 

2.2.3. 1  Implementation  of  the  DPDE  Pair  Style 

Implementation  of  constant  energy  DPD  with  the  VV  integration  scheme  requires  a  new  pair 
style  to  be  defined,  which  computes  the  conservative,  dissipative,  and  random  forces  as  well  as 
the  change  in  particle  internal  conductive  energy  and  mechanical  energy  from  one  integration 
step  to  the  next.  The  new  PairDPDE  class  (LAMMPS  source  files:  pair  dpde.h  and 
pair_dpde.cpp)  is  specified  through  the  LAMMPS  input  file  via  the  pair_style  command 

pair_style  dpde  <kappa  flag>  <cutoff>  <random  number  seed> 
pair_coeff  i  j  ADPD  sigmaij  kappa  <cutoff> 

An  energy-independent  or  energy-dependent  (10)  kappa  model  is  specified  in  the  pair_style 
command  by  setting  the  kappaflag  to  0  or  1 ,  respectively.  In  the  energy-independent  model 
(kappaflag  =  0),  Ky  is  explicitly  given  as  a  pair  coefficient  in  the  LAMMPS  input  file.  In  the 
energy-dependent  model  (kappa_flag  =  1),  ic,,  is  given  by  the  equation 
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(6) 


KU  = 


Kq 

fcB 


2 


where  the  Ko  pair  coefficient  is  specified  in  the  LAMMPS  input  file,  U;  is  the  total  internal  energy 
of  particle  i,  and  kB  is  Boltzmann’s  constant. 

The  PairDPDE  pair  style  class  derives  from  the  Pair  parent  class  and  is  very  similar  to  and 
modeled  after  the  existing  PairDPD  pair  style  class,  except  that  it  contains  additional  per-atom 
arrays  that  update  the  DPD  particle  internal  temperature  (dpdTheta),  conductive  energy  (uCond), 
mechanical  energy  (uMech),  heat  capacity  (cv),  and  density  (rho),  as  well  as  the  differences  in 
the  DPD  particle’s  conductive  energy  (duCond)  and  mechanical  energy  (duMech)  between  two 
subsequent  integration  steps.  The  main  differences  between  the  native  PairDPD  and  modified 
PairDPDE  classes  are  summarized  in  appendix  C. 

2. 2. 3. 2  Implementation  of  DPDE  Fix  Commands 

Implementation  of  constant  energy  DPD  with  the  VV  integration  scheme  also  requires  a  new  fix 
class  to  be  defined,  which  accounts  for  the  integration  of  the  conductive  and  mechanical  energy 
and  applying  a  mesoparticle  equation  of  state.  The  new  fix  is  defined  as  the  FixDPDE  class 
(LAMMPS  source  files:  fix_dpde.h  and  fix_dpde.cpp)  and  is  specified  in  the  LAMMPS  input  file 
via  the  fix  command 

fix  1  all  dpde 

The  FixDPDE  is  very  similar  to  and  modeled  after  the  existing  FixNVE  class,  except  that  it 
contains  the  integration  of  the  conductive  and  mechanical  energy.  The  main  differences  between 
the  native  FixNVE  and  modified  FixDPDE  classes  are  summarized  in  appendix  D. 

Once  the  updated  conductive  and  mechanical  energy  is  computed,  the  DPD  particle  internal 
temperatures  are  updated  through  the  mesoparticle  equation  of  state  that  relates  the  total  particle 
internal  particle  energy  to  the  internal  temperature.  The  mesoparticle  equation  of  state  must  be 
specified  with  a  separate  fix  command.  Currently,  only  one  choice  of  the  mesoparticle  equation 
of  state  has  been  implemented.  It  relates  the  total  internal  energy  to  the  internal  temperature 
through  the  heat  capacity  according  to  the  relation 

ui  =  Cv,i@i  (7) 

The  new  fix  is  defined  as  the  FixEOScv  class  (LAMMPS  source  files:  fix_eos_cv.h  and 
fix_eos_cv.cpp )  and  is  specified  in  the  LAMMPS  input  file  via  the  fix  command: 

fix  1  all  eos/cv  <cv> 

where  cv  is  the  value  of  the  heat  capacity.  The  FixEOScv  class  is  summarized  in  appendix  E. 
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2.2.4  Implementation  of  DPD-E  With  the  VV-SSA  Integration  Scheme 

Implementation  of  the  isoenergetic  DPD  using  the  VV-SSA  integration  scheme  requires  a 
splitting  of  the  stochastic  and  deterministic  EOMs.  The  conservative  force  is  computed  within 
the  dpde/conservative  pair  style  and  is  integrated  deterministically  through  the  VV  integration 
scheme  with  the  fix  dpde/shardlow  command.  The  random  and  dissipative  forces  are  computed 
and  integrated  through  the  SSA  within  the  fix  dpde/shardlow  prior  to  the  deterministic 
integration  of  the  conservative  force.  The  new  pair  style  and  fix  commands  are  presented 
together  since  the  classes  interact  with  one  another  to  perfonn  all  stages  of  the  VV-SSA 
integration.  Additional  implementation  details  of  the  SSA  can  be  found  elsewhere  (9). 

2.2.4. 1  Implementation  of  the  DPDE/Conservative  Pair  Style 

Implementation  of  constant  energy  DPD  using  the  VV-SSA  integration  scheme  requires  a  new 
pair  style  that  computes  the  conservative  force  to  be  defined.  The  new  PairDPDEConservative 
class  (LAMMPS  source  files:  pa ir  dpde  cons erva ti ve. h  and  pairdpdeconservative. cpp )  is 
specified  in  the  LAMMPS  input  file  via  the  pair  style  command 

pair  style  dpde/conservative  <kappa  flag>  <cutoff>  <random  number  seed> 
pair_coeff  i  j  ADPD  sigmaij  kappa  <cutoff> 

An  energy-independent  or  energy-dependent  (10)  kappa  model  is  specified  in  the  pair  style 
command  by  setting  the  kappa  flag  to  0  or  1 ,  respectively.  In  the  energy-independent  model 
(kappa_flag  =  0),  iq,  is  explicitly  given  as  a  pair  coefficient  in  the  LAMMPS  input  file.  In  the 
energy-dependent  model  (kappa_flag  =  1),  iq,  is  related  to  the  k0  pair  coefficient  as  given  in 
equation  6. 

The  PairDPDEConservative  pair  style  class  derives  from  the  Pair  parent  class  and  is  very  similar 
to  and  modeled  after  the  existing  PairDPD  pair  style  class.  The  main  difference  is  that  random 
and  dissipative  computations  are  removed,  while  the  required  data  structures  to  compute  the 
DPD  particle  attributes  are  included.  The  main  differences  between  the  PairDPD  and 
PairDPDEConservative  classes  are  summarized  in  section  2. 2.4. 2  and  are  explicitly  shown  in 
appendix  F. 

2.2. 4.2  Implementation  of  DPDE/Shardlow  Fix  Command 

Implementation  of  constant  energy  DPD  using  the  VV-SSA  integration  also  requires  a  new  fix 
class  to  be  defined,  which  accounts  for  the  stochastic  integration  of  the  random  and  dissipative 
forces,  the  deterministic  integration  of  the  conservative  force,  and  the  integration  of  the 
conductive  and  mechanical  energy  for  constant  energy  simulations.  The  new  fix  is  defined  as  the 
FixDPDEShardlow  class  (LAMMPS  source  files:  pair_dpde_shardlow.h  and pair_dpde_ 
shardlow.cpp )  and  is  specified  in  the  LAMMPS  input  file  via  the  fix  command: 

fix  1  all  dpde/shardlow 
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The  FixDPDEShardlow  class  is  very  similar  to  the  existing  FixNVE  class,  except  that  it  contains 
an  additional  function  to  integrate  the  stochastic  equations  of  motion,  the  conductive  energy,  and 
the  mechanical  energy.  These  integrations  are  handled  via  a  function  called  stochastic  integrate 
within  the  FixDPDEShardlow  class.  The  DPD  particle  internal  temperatures  are  updated  through 
the  mesoparticle  equation  of  state  that  relates  the  total  particle  internal  particle  energy  to  the 
internal  temperature.  The  mesoparticle  equation  of  state  must  be  specified  with  a  separate  fix 
command.  Currently,  only  one  choice  of  the  mesoparticle  equation  of  state  has  been 
implemented  and  can  be  specified  in  the  LAMMPS  input  file  via  the  fix  command 

fix  1  all  eos/cv  <cv> 

The  main  differences  between  the  native  FixNVE  and  modified  FixDPDEShardlow  classes  are 
summarized  in  appendix  G. 


3.  Example  Input  Files  for  VV  and  V V-SSA  Integration 


Example  input  files  running  a  DPD-E  simulation  in  LAMMPS  with  the  VV  and  VV-SSA 
integration  schemes  are  provided  in  figures  1  and  2,  respectively. 
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Figure  1.  Example  input  file  for  simulating  a  DPD  fluid  under  isoenergetic  conditions  with  the  VV  integration 
scheme. 

#  Input  File  for  DPD  fluid  under  isoenergetic  conditions  using  the  VV  integration  scheme 
boundary  p  p  p 

units  metal  #ev,  ps 

atom_style  dpd 

read_data  initial. conf.DPDfluid 

communicate  single  vel  yes 

mass  1  125.9 

pair_style  dpde  1  8.60  234324 

pair_coeff  1  1  0.0752  0.0223  4.55E-05  8.60 

neighbor  2.0  bin 

neighjmodify  every  1  delay  0  check  no  once  no 

#  Time  in  ps  for  metal  units 

timestep  0.001 

compute  dpdU  all  dpd 

variable  totEnergy  equal  pe+ke+c_dpdU[3] 

thermo  1 

thermo_style  custom  step  temp  pe  ke  c_dpdU[l]  c_dpdU[2]  c_dpdU[3]  v_totEnergy  c_dpdll[4] 
thermo_modify  format  float  %15.10E 

fix  1  all  dpde 

fix  2  all  eos/cv  0.00517041 

run  100 
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Figure  2.  Example  input  file  for  simulating  a  DPD  fluid  under  isoenergetic  conditions  with  the  VV-SSA 
integration  scheme. 


#  Input  File  for  DPD  fluid  under  isoenergetic  conditions  using  the  VV-SSA  integration  scheme 
boundary  p  p  p 


units 

metal  #  ev,  ps 

atom_style 

dpd 

read_data 

initial.  conf.DPDfluid 

communicate  single  vel  yes 

mass 

1  125.9 

pairstyle 

dpde/conservative  1  8.60  234324 

pair_coeff 

1  1  0.0752  0.0223  4.55E-05  8.60 

neighbor 

2.0  bin 

neigh_modify  every  1  delay  0  check  no  once  r 

#  Time  in  ps  for  metal  units 

timestep 

0.001 

compute 

dpdU  all  dpd 

variable 

totEnergy  equal  pe+ke+c_dpdU[3] 

thermo 

1 

thermo_style  custom  step  temp  pe  ke  c_dpdll[l]  c_dpdU[2]  c_dpdU[3]  v_totEnergy  c_dpdll[4] 
thermo_modify  format  float  %15.10E 

fix  1  all  dpde/shardlow 

fix  2  all  eos/cv  0.00517041 

run  100 


The  example  input  files  contain  the  commands  to  simulate  a  10,125  particle  system  in  a  cubic 
box  with  volume  129.0  A3  at  a  temperature  of  300  K  under  isoenergetic  conditions  with  a  time 
step  of  0.001  ps.  The  eos/cv  equation  of  state  is  specified  along  with  an  energy-dependent  kappa 
model.  The  neighbor  lists  are  reconstructed  after  every  time  step,  and  the  thermodynamics  are 
output  every  0.001  ps.  The  total  particle  internal  energies  and  internal  temperature  are  included 
in  the  thennodynamics  via  the  compute  dpd  command.  The  DPD-E  parameters  are  summarized 
in  table  1 . 
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Table  1.  Model  and  system  parameters  for  DPD-E  simulations 
of  the  DPD  fluid  model. 


System  Property 

Real  Units 

Pdpd 

4.72  x  10'3  A3 

Tsp  —  Tinit 

300  K 

mass 

125.9  g  mol1 

Adpd 

7.52  x  10'2  eV  A'1 

T'cut 

8.6  A 

o 

2.23  x  10'2  eV  ps1/2  A"1 

K0 

4.55  x  10'5  ps"1 

cv 

5.17  x  10'3  eVK1 

4.  Results 


The  example  benchmarks  described  in  section  5  are  compared  in  terms  of  performance  and 
accuracy  in  Lisal  et  al.  (5),  where  DPD  fluid  calculations  were  perfonned  under  isoenergetic 
conditions  with  10,125  particles  at  time  steps  ranging  from  0.001  to  0.4  ps  for  a  total  of  1  ns  of 
simulation  time.  Lisal  et  al.  (5)  shows  that  the  VV-SSA  integrator  is  considerably  more  stable 
than  the  VV  integrator,  where  comparable  accuracy  is  achieved  with  time  steps  that  are  10-100 
times  longer  than  the  traditional  VV  integrator.  However,  the  VV-SSA  implementation  is  shown 
to  be  more  computationally  expensive  than  the  serial  VV  implementation  on  a  per  time  step  basis 
and  carries  additional  parallelization  overhead.  Analysis  of  the  trade-off  between  computation 
performance  and  time-step  size  stability  showed  that  the  VV-SSA  integration  scheme  can 
decrease  the  overall  time-to-solution  by  a  factor  of  10-100  for  a  DPD-E  simulation,  thereby 
justifying  practical  and  regular  implementation  of  the  approach. 


5.  Conclusions 


In  this  report  we  have  provided  the  user  with  the  following  infonnation  and  resources: 

•  A  detailed  description  of  the  implementation  of  the  DPD-E  method  using  the  VV  and 
VV-SSA  integration  schemes  into  LAMMPS. 

•  Documentation  that  highlights  the  major  modifications  to  the  source  code. 

•  Example  input  scripts  to  run  the  software. 

•  A  benchmark  study  to  give  an  indication  of  the  performance  and  accuracy  of  DPD-E  with 
the  VV  and  VV-SSA  integration  schemes. 
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Appendix  A.  Code  Differences  Between  the  Atom  Class  Files  as  Compared  to 

the  Native  LAMMPS  Code 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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diff  --git  a/src/atom.h  b/src/atom.h 
index  cc044f1..1d122ec  100644 
—  a/src/atom.h 
+++  b/src/atom.h 

@@  -50,6  +50,10  @@  class  Atom  :  protected  Pointers  { 
int  *type,*mask; 
imageint  *image; 
double  **x,**v,**f; 

+  double  *uCond,  *uMech; 

+  double  *duCond,  *duMech; 

+  double  *dpdTheta; 

+  int*eos; 

tagint  *molecule; 
int  *molindex,*molatom; 


diff  --git  a/src/atom.cpp  b/src/atom.cpp 
index  86ae09b..888aebe  100644 
—  a/src/atom.cpp 
+++  b/src/atom.cpp 

@@-72,6  +72,10  @@  Atom  "Atom  (LAMM  PS  *lmp)  :  Pointers(lmp) 
type  =  mask  =  NULL; 
image  =  NULL; 
x  =  v  =  f  =  NULL; 

+  uCond  =  uMech  =  NULL; 

+  duCond  =  duMech  =  NULL; 

+  dpdTheta  =  NULL; 

+  eos  =  NULL; 

molecule  =  NULL; 
molindex  =  molatom  =  NULL; 
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Appendix  B.  Summary  of  LAMMPS  Code  Differences  Between  the 
AtomVecDPD  Class  as  Compared  to  the  Native  AtomVecAtomic  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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AtomVecDPD  Constructor: 

O  Initialized  uCond,  uMech,  duCond,  duMech  and  eos  pointers  to  NULL 

O  Updated  comm_x_only  and  comm._f_only  to  1  to  allow  for  forward  communication  of  rho,  cv  and  dpdTheta 
O  Updated  data  buffer  size  definitions  for  the  forward  and  border  communications  (size_forward=9,  size_border=12)  to 
account  for  cv,  dpdTheta,  eos,  uCond  and  uMech 

O  Updated  the  number  of  columns  in  the  data  file  (size_data_atom  =  6).  The  data  file  contains  the  following  six  columns: 

■  Column  1:  particle  id 

■  Column  2:  particle  type 

■  Column  3:  Internal  temperature  (0.)  of  particle  i 

■  Columns  4-6:  Cartesian  coordinates  (x,  y,  z)  of  particle  i 
O  Updated  the  column  the  x  data  is  located  (xcol_data  =  4) 

O  Set  the  rho_flag  to  be  on  (atom->rho_flag  =  1) 
void  AtomVecDPD::grow(int  n) 

O  Added  memory  allocation  for  rho,  drho,  cv,  dpdTheta,  eos,  uCond,  uMech,  duCond  and  duMech 
void  AtomVecDPD::grow_reset() 

O  Added  pointers  for  rho,  drho,  cv,  dpdTheta,  eos,  uCond,  uMech,  duCond  and  duMech 
void  AtomVecDPD::copy(int  i,  int  j,  int  delflag) 

O  Copies  particle  i  rho,  drho,  cv,  dpdTheta,  eos,  uCond,  uMech  to  particle  j 
int  AtomVecDPD: :pack_comm(int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  forward  communication  buffer 
int  AtomVecDPD::pack_comm_vel(int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  forward  communication  buffer,  along  with  velocities 
void  AtomVecDPD::unpack_comm(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  forward  communication  buffer 
void  AtomVecDPD::unpack_comm_vel(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  forward  communication  buffer,  along  with  velocities 
int  AtomVecDPD::pack_border(int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  border  communication  buffer, 
int  AtomVecDPD: :pack_border_vel (int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  border  communication  buffer,  along  with  velocities 
void  AtomVecDPD::pack_comm_hybrid(int  n,  int  *list,  double  *buf) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  for  the  hybrid  atom  style  communication  buffer 
void  AtomVecDPD::pack_border_hybrid(int  n,  int  *list,  double  *buf) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  hybrid  atom  style  border  communication  buffer 
void  AtomVecDPD::unpack_border(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  the  border  communication  buffer 
void  AtomVecDPD::unpack_border_vel(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  the  border  communication  buffer,  along  with  velocities 
void  AtomVecDPD::unpack_comm_hybrid(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  the  hybrid  atom  style  communication  buffer 
void  AtomVecDPD::unpack_border_hybrid(int  n,  int  first,  double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  the  hybrid  atom  style  border  communication  buffer 
int  AtomVecDPD: :pack_exchange(int  i,  double  *buf) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  exchange  communication  buffer, 
int  AtomVecDPD: :unpack_exchange(double  *buf) 
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O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  exchange  communication  buffer, 
int  AtomVecDPD::size_restart() 

O  Adjusted  the  size  of  the  restart  data  to  include  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech 
■  Changed  from  1 1  to  1 7 
int  AtomVecDPD::pack_restart(int  i,  double  *buf) 

O  Packs  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  in  the  restart  buffer 
int  AtomVecDPD::unpack_restart(double  *buf) 

O  Unpacks  rho,  cv,  dpdTheta,  eos,  uCond  and  uMech  from  the  restart  buffer 
void  AtomVecDPD::create_atom(int  itype,  double  *coord) 

O  Creates  a  new  particle  and  sets  default  rho,  drho,  cv,  dpdTheta,  eos,  uCond,  uMech,  duCond  and  duMech  values 
void  AtomVecDPD::data_atom(double  *coord,  tagint  imagetmp,  char  **values) 

O  Read  internal  temperature  (dpdTheta)  from  data  file;  initialize  rho,  cv,  uCond[i]  and  uMech[i]  to  zero  and  eos  to  -1 
int  AtomVecDPD::data_atom_hybrid(int  nlocal,  char  **values) 

O  Read  internal  temperature  (dpdTheta)  from  hydrid  atom  style,  data  file 
void  AtomVecDPD::pack_data(double  **buf) 

O  Packs  internal  temperature  (dpdTheta)  into  the  data  file  buffer 
int  AtomVecDPD::pack_data_hybrid(int  i,  double  *buf) 

O  Packs  internal  temperature  (dpdTheta)  into  the  hydrid  atom  style,  data  file  buffer 
void  AtomVecDPD::write_data(FILE  *fp,  int  n,  double  **buf) 

O  Writes  internal  temperature  (dpdTheta)  to  the  data  file 
int  AtomVecDPD::write_data_hybrid(FILE  *fp,  double  *buf) 

O  Writes  internal  temperature  (dpdTheta)  to  the  hybrid  atom  info  to  data  file 
bigint  AtomVecDPD::memory_usage() 

O  Added  memory  checks  for  rho,  drho,  cv,  dpdTheta,  eos,  uCond,  uMech,  duCond  and  duMech 
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Appendix  C.  Summary  of  LAMMPS  Code  Differences  Between  the  PairDPDE 
Class  as  Compared  to  the  Native  PairDPD  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 


•  PairDPDE  constructor:  PairDPDE::PairDPDE(LAMMPS  *lmp)  :  Pair(lmp) 

O  Set  the  size  of  the  pair  forward  communication  to  3  (duCond^  duMech,,  dpdTheta^ 

O  Set  the  size  of  the  pair  reverse  communication  to  2  (duCond,,  duMech,) 


PairDPDE  destructor:  PairDPDE::~PairDPDE() 

O  Deallocate  kappa  array;  remove  gamma  deallocation 
void  PairDPDE::compute(int  eflag,  int  vflag) 

O  Initialize  duCond  and  duMech  by  setting  the  arrays  to  0.0  followed  by  a  forward  communication 

O  Define  wR  and  wD:  0)D (r)  =  [ojfl  (r)]2  =  j(^  rc)  •  r  <  rc 

(  0  ,  r  >rc 

where  rc  is  the  cutoff  radius  for  the  pair  of  interacting  particles. 

NOTE:  In  PairDPD  class,  wd  definition  differs 

1 (l  .  l\ 

O  Compute  thetajj:  t —  —  I  — — h  — —  I 

J  2  \di  OjJ 


O  Compute  gammajj:  Yij  —  —  — 


ij 


O  Compute  the  conservative,  dissipative  and  random  forces 

F  =  Fc  +  F0  +  Fs 
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where  is  a  random  number  sampling  a  Gaussian  distribution  and  At  is  the  timestep  used  in  the  simulation 
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O  Compute  kappa_ij:  Ky  or  Kij  —  —  ^ — ~ — J 
O  Compute  alpha_ij:  CCij  = 

(  1  1\ 

O  Compute  mu_ij:  j U//  —  I - 

J  \m  ™-j) 

O  Compute  duMech  and  duMech 
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O  Compute  duCondi  and  duCond 
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O  Compute  the  interaction  energy  between  the  particle  pair 

■  evdwl  =  0.5*a0[itype][jtype]*cut[itype][jtype]  *  wd; 
O  Reverse  communication  of  duCond  and  duMech 
void  PairDPDE::allocate() 

O  Allocate  memory  for  kappa  array;  Remove  gamma  allocation 
void  PairDPDE::settings(int  narg,  char  **arg) 
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O  Read  3  arguments  associated  with  pair  style  (kappa_flag,  cut_global  and  seed);  Remove  temperature 
O  Ensure  kappa_flag  is  set  to  either  0  or  1 

•  void  PairDPDE::coeff(int  narg,  char  **arg) 

O  Reads  in  3-4  three  pair  coefficients 

■  A 

'tipd 

■  Sigma  (0£y);  removed  gamma 

■  Kappajj  ( Kij )  or  KappaO  (Kq) 

■  cutoff  (optional) 

•  void  PairDPDE::init_style() 

O  Ensures  that  the  Newton  pair  is  set  to  “on"  for  DPD  calculations 

•  double  PairDPDE::init_one(int  i,  int  j) 

O  Initializes  the  pair  coefficients  for  all  particle  pairs;  Add  kappa,  sigma  and  remove  gamma 

•  void  PairDPDE::write_restart(FILE  *fp) 

O  Updates  the  pair  coefficient  information  to  the  restart  file;  Add  kappa,  sigma  and  remove  gamma 

•  void  PairDPDE::read_restart(FILE  *fp) 

O  Read  pair  coefficients  from  restart  file  and  broadcast  across  processors;  Add  kappa,  sigma  and  remove  gamma 

•  void  PairDPDE::write_restart_settings(FILE  *fp) 

O  Removed  the  temperature  variable 

•  void  PairDPDE::read_restart_settings(FILE  *fp) 

O  Removed  temperature  variable 

•  void  PairDPD::write_data(FILE  *fp) 

O  Removed  gamma;  Added  kappa  and  sigma  to  list  of  variable  to  write  to  the  data  file 

•  void  PairDPD::write_data_all(FILE  *fp) 

O  Removed  gamma;  Added  kappa  and  sigma  to  list  of  variable  to  write  to  the  data  file 

•  double  PairDPDE::single(int  i,  int  j,  int  itype,  int  jtype,  double  rsq,  double  factor_coul,  double  factor_dpd,  double  &fforce) 

O  Updated  definition  of  wr  and  wd 

•  int  PairDPDE::pack_comm(int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  the  buffer  for  forward  communication  (duCond,  duMech  and  dpdTheta) 

•  void  PairDPDE::unpack_comm(int  n,  int  first,  double  *buf) 

O  Unpacks  the  buffer  for  forward  communication  (duCond,  duMech  and  dpdTheta) 

•  int  PairDPDE::pack_reverse_comm(int  n,  int  first,  double  *buf) 

O  Packs  the  buffer  for  reverse  communication  (duCond  and  duMech) 

•  void  PairDPDE::unpack_reverse_comm(int  n,  int  *list,  double  *buf) 

O  Unpacks  the  buffer  for  reverse  communication  (duCond  and  duMech) 
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Appendix  D.  Summary  of  LAMMPS  Code  Differences  Between  the  FixDPDE 
Class  as  Compared  to  the  Native  FixNVE  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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FixDPDE  constructor 


O  Ensure  3  keywords  are  given  in  the  fix  dpde  command 

•  int  FixDPDE::setmask() 

O  Removed  all  reference  to  RESPA 

•  void  FixDPDE::init() 

O  Ensure  that  an  equation  of  state  is  specified  in  another  fix;  otherwise,  print  error  message  and  stop 
O  Remove  all  references  to  RESPA 

•  void  FixDPDE::initial_integrate(int  vflag) 

O  Integrate  uCond  and  uMech  by  the  relations: 

uCondi  =  uCondi  +  14  *  dt  *  duCond; 
uMechi  =  uMechi  +  14  *  dt  *  duMechi 

•  void  FixDPDE::final_integrate() 

O  Integrated  uCond  and  uMech  by  the  relations: 

uCondi  =  uCondi  +  14  *  dt  *  duCondi 
uMechi  =  uMechi  +  14  *  dt  *  duMechi 
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Appendix  E.  Summary  of  the  FixEOScv  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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FixEOScv  constructor 


O  Ensure  4  keywords  are  specified  in  the  fix  eos/cv  command 
O  Read  the  heat  capacity,  cv 

•  int  FixEOScv::setmask() 

O  Set  the  mask  to  POSTJ NTEGRATE  and  END_OF_STEP 

•  void  FixEOScv::init() 

O  Initialize  the  cv  array 
O  Initialize  the  eos  array 

O  Initialize  uCond  and  uMech  to  be  *  cv[i]  *  dpdTheta[i] 

•  void  FixEOScv: :post_integrate() 

O  Apply  the  EOS  to  compute  the  DPD  particle  internal  temperature.  The  eos/cv  defines  the  following  relation: 
0i  =  (uCondi  +  uMech,)  /  Cvi 
where  0j  is  the  internal  particle  temperature  (dpdTheta). 

•  void  FixEOScv::end_of_step() 

O  Apply  the  EOS  to  compute  the  DPD  particle  internal  temperature.  The  eos/cv  defines  the  following  relation: 
0i  =  (uCondi  +  uMech,)  /  Cvi 

where  0,  is  the  internal  particle  temperature  (dpdTheta). 
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Appendix  F.  Summary  of  LAMMPS  Code  Differences  Between  the  PairDPDE 
Conservative  Class  as  Compared  to  the  Native  PairDPD  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 


29 


PairDPDEConservative  destructor:  PairDPDE::~PairDPDE() 

O  Deallocate  kappa  array;  remove  gamma  deallocation 


void  PairDPDEConservative::compute(int  eflag,  int  vflag) 

O  Removed  dissipative  and  random  force  calculations 


O  Define  wR  and  wD:  (jl)D  (r)  —  [(x)^(r)]2 


r  <rc 
r  >rc 


where  rc  is  the  cutoff  radius  for  the  pair  of  interacting  particles. 
NOTE:  In  PairDPD  class,  wd  definition  differs 
O  Compute  the  conservative  force:  Ffj  —  ADpDaR(r) 

O  Compute  the  interaction  energy  between  the  particle  pair 

■  evdwl  =  0.5*a0[itype][jtype]*cut[itype][jtype]  *  wd; 


void  PairDPDEConservative: :allocate() 

O  Allocate  memory  for  kappa  array;  Remove  gamma  allocation 


void  PairDPDEConservative::settings(int  narg,  char  **arg) 

O  Read  3  pair  style  arguments  (kappa_flag,  cut_global  and  seed);  Remove  the  temperature  argument 
O  Ensure  kappa_flag  is  set  to  either  0  or  1 


void  PairDPDEConservative::coeff(int  narg,  char  **arg) 

O  Reads  in  3-4  three  pair  coefficients 

■  A 

^DPD 

■  Sigma  (0£y);  removed  gamma 

■  Kappajj  {K(j)  or  KappaO  (JCq) 

■  cutoff  (optional) 
void  PairDPDEConservative: :init_style() 

O  Ensures  that  ghost  velocities  are  stored  and  Newton  pair  is  set  to  “on"  for  DPD  calculations 
double  PairDPDEConservative::init_one(int  i,  int  j) 

O  Initializes  the  pair  coefficients  for  all  particle  pairs;  Add  kappa,  sigma  and  remove  gamma 
void  PairDPDEConservative::write_restart(FILE  *fp) 

O  Updates  the  pair  coefficient  information  to  the  restart  file;  Add  kappa,  sigma  and  remove  gamma 
void  PairDPDEConservative::read_restart(FILE  *fp) 

O  Reads  pair  coefficients  from  restart  file  and  broadcast  across  processors;  Add  kappa,  sigma  and  remove  gamma 
void  PairDPDEConservative::write_restart_settings(FILE  *fp) 

O  Remove  the  temperature  variable 
void  PairDPDEConservative::read_restart_settings(FILE  *fp) 

O  Remove  temperature  variable 
void  PairDPDEConservative::write_data(FILE  *fp) 

O  Removed  gamma;  Added  kappa  and  sigma  to  list  of  variable  to  write  to  the  data  file 
void  PairDPDEConservative: :write_data_all(FILE  *fp) 

O  Removed  gamma;  Added  kappa  and  sigma  to  list  of  variable  to  write  to  the  data  file 


double  PairDPDEConservative::single(int  i,  int  j,  int  itype,  int  jtype,  double  rsq,  double  factor_coul,  double  factor_dpd,  double  &fforce) 
O  changed  definition  of  wd  and  wr  to  be  consistent  with  manuscript 
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Appendix  G.  Summary  of  LAMMPS  Code  Differences  Between  the 
FixDPDEShardlow  Class  as  Compared  to  the  Native  FixNVE  Class 


This  appendix  appears  in  its  original  form,  without  editorial  change. 


•  FixDPDEShardlow  Constructor 

O  Ensure  3  keywords  are  given  in  the  fix  dpde/shardlow  command 

O  Set  the  size  of  the  fix  forward  communication  to  1 0  (dvSSA[i][0-2],  v[i] [0-2] ,  duCondp],  duMech[i],  uCond[i],  uMech[i]) 

O  Set  the  size  of  the  fix  reverse  communication  to  5  (dvSSA[i][0-2], duCondp],  duMech[i]) 


•  int  FixDPDEShardlow::setmask() 

O  Remove  all  references  to  RESPA 

•  void  FixDPDEShardlow::init() 

O  Ensure  that  an  equation  of  state  is  specified  in  another  fix;  otherwise,  print  error  message  and  stop 
O  Set  the  pairDPDE  pointer  to  access  the  pairDPDEConservative  class  data  structures 

•  void  FixDPDEShardlow::initial_integrate(int  vflag) 

O  Added  a  shardlow_integrate()  function  prior  to  the  velocity  integration 
O  Integrate  the  velocity  and  position  at  a  half  time  step  (Same  as  FixNVE) 

•  void  FixDPDEShardlow: :final_integrate() 

O  Integrate  the  velocities  at  a  half  time  step  (Same  as  FixNVE) 

•  void  FixDPDEShardlow::shardlow_integrate():  Member  function  to  perform  the  stochastic  integration  of  velocities 

O  Compute  the  length  of  the  bounding  box  dimensions  (bbx,  bby,  bbz) 

■  Ensure  all  bounding  box  dimensions  are  larger  than  interaction  cutoff  radius;  Print  error  if  not  satisfied 
O  Allocate  memory  for  velocity  (dvSSA[j][0-2]) 

O  Initialize  (dvSSA[j][0-2],  duCondp],  duMechp])  arrays  and  forward  communicate  to  ensure  all  processors  are  initialized 
O  Allocate  memory  for  the  particle  pair  active  interaction  regions  (jbin) 

■  Count  the  number  of  pairs  (count),  then  allocate  active  interaction  region  bins,  jbin[count] 

■  Assign  each  pair  to  an  active  interaction  region  depending  on  the  coordinates  of  the  pair.  One  atom  will  be  in 
central  box;  assign  the  AIR  based  on  the  2nd  atom  .which  may  be  in  the  central  box  or  a  neighboring  box 

O  Loop  over  the  eight  active  interaction  regions  (Stages) 

■  If  the  pair  lies  in  the  current  active  interaction  region,  then  . . . 

•  Compute  the  pair  separation  distance,  rr  and  check  if  it  is  less  than  the  interaction  cutoff 
O  Store  the  current  velocities  (vxOi,  vyOi,  vzOi,  vxOj,  vyOj,  vzOj) 


O 

O 

O 

O 

O 

O 

O 

O 


P,  P/ 

Compute  the  velocity  difference  between  i  and  j:  V  ••  ^ - =  vi "  vi 

mi  rrij 


Compute  the  dot  product,  dot  =  ry  •  vy 
Compute  wR  and  wD:  (jl)D  (r)  —  [(x)^(r)]2  — 


(i-^)  ,r<rc 
l  0  ,  r  >rc 

where  rc  is  the  cutoff  radius  for  the  pair  of  interacting  particles 
Apply  EOS  to  compute  the  current  internal  temperature. 


Compute  thetajj:  I 
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Obtain  a  random  number  from  a  Gaussian  distribution 
Compute  the  momentum  (velocity)  update 
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O  Re-Compute  the  velocity  difference  between  i  and  j:  V  ••  ^ - =  vi  ■  vj 
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O  Compute  the  momentum  (velocity)  update 


P,  <-  P, 


At 


YiP 


^  \  +  —yjj(oD A t 
2  7 


( 


ru 

—  ■  V.. 

r  lJ 

\r‘i  J 


rJL  +  Z.aS(,rJLfi; 

r..  2  r... 


+  o-,®  c, 


Va7 


^  2 


P^P  + 


Ar 


1  +  —  Yu®0  At 
2  7 


ILV 

ii 

r. 

v  y 


y 


ry  2  'i/ 


r//  Va7 

'<¥»  ?u  0 

^  2 


O  Compute  a  new  random  number, 

K0  / Uj+UA 2 

O  Compute  kappajj:  Kij  —  —  ^ — - — J  (for  T-dependent  model) 

O  Compute  alpha_ij:  CCy  =  ^Zk^Kjj 
O  Compute  the  conductive  energy  update 


„  cond  .  ,  cond  . 
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O  Compute  the  mechanical  energy  update 
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•  Reverse  communicate  the  velocity  updates. 

•  Update  the  velocities  and  energies  (v[j][0-2],  uCondp],  uMech[j]) 

•  Reinitialize  the  arrays  to  zero  (dvSSA[j][0-2],  duCond[j],  duMechp]) 

•  Forward  communicate  the  velocity/energy  updates 

■  Delete  the  memory  allocations  for  jbin  and  dvSSA 
int  FixDPDEShardlow::sort_bin(double  rx,  double  ry,  double  rz) 

O  Compare  the  particle  position  to  the  reference  points  of  the  bounding  box 

■  Determine  whether  the  point  lies  inside  the  bounding  box,  or  outside  the  bounding  box 

■  Return  a  number  to  indicate  the  active  interaction  region 

•  0  =  Inside  bounding  box 

•  1  =  Top 

•  2  =  Right 
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•  3  =  Top-right  and  Bottom-right 

•  4  =  Back 

•  5  =  Top-back  and  Bottom-back 

•  6  =  Left-back  and  Right-back 

•  7  =  Back  Corners 

•  int  FixDPDEShardlow::pack_comm(int  n,  int  *list,  double  *buf,  int  pbc_flag,  int  *pbc) 

O  Packs  the  buffer  for  forward  communication  (dvSSA,  v,  duCond,  duMech,  uCond,  uMech) 

•  void  FixDPDEShardlow::unpack_comm(int  n,  int  first,  double  *buf) 

O  Unpacks  the  buffer  for  forward  communication  (dvSSA,  v,  duCond,  duMech,  uCond,  uMech) 

•  int  FixDPDEShardlow::pack_reverse_comm(int  n,  int  first,  double  *buf) 

O  Packs  the  buffer  for  reverse  communication  (dvSSA,  duCond,  duMech) 

•  void  FixDPDEShardlow::unpack_reverse_comm(int  n,  int  *list,  double  *buf) 

O  Unpacks  the  buffer  for  reverse  communication  (dvSSA,  duCond,  duMech) 
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